library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(ggplot2)
library(sdmTMB)
library(tidyr)

Mesh

# 300 knots
plot(mesh)

Model specifications

m1 <- sdmTMB(
  bio.adult ~ 0 + YearF + GearCat + s(Month, bs = "cc", k = 10),
  data = dat,
  mesh = mesh,
  offset = log(dat$SweptArea),
  family = tweedie(link = "log"),
  spatial = "on",
  time = "Year",
  spatiotemporal = "IID"
)
## Warning in checkDepPackageVersion(dep_pkg = "TMB"): Package version inconsistency detected.
## glmmTMB was built with TMB version 1.8.1
## Current TMB version is 1.9.0
## Please re-install glmmTMB from source or restore original 'TMB' package (see '?reinstalling' for more information)
## Warning: The model may not have converged: non-positive-definite Hessian matrix.

Some diagnostics

dat$resids <- residuals(m1) # randomized quantile residuals

hist(dat$resids)

qqnorm(dat$resids[is.finite(dat$resids)])
qqline(dat$resids[is.finite(dat$resids)])

# check spatial autocorrelation
ggplot(filter(dat,Year%%5==1), aes(lon, lat, col = resids)) +
  scale_colour_gradient2() +
  geom_point() + geom_sf(data=map,inherit.aes = FALSE)+xlim(-12,30)+ylim(47,62)+
  facet_wrap(~Year,ncol=2) 
## Warning: Removed 758 rows containing missing values (geom_point).

Terms

m1
## Spatiotemporal model fit by ML ['sdmTMB']
## Formula: bio.adult ~ 0 + YearF + GearCat + s(Month, bs = "cc", k = 10)
## Mesh: mesh
## Time column: Year
## Data: dat
## Family: tweedie(link = 'log')
##  
##               coef.est coef.se
## YearF1967         0.00     NaN
## YearF1968         0.00     NaN
## YearF1970        11.34     NaN
## YearF1971        10.15     NaN
## YearF1972        10.10     NaN
## YearF1973        10.01     NaN
## YearF1974         9.99     NaN
## YearF1975         9.08     NaN
## YearF1976        10.04     NaN
## YearF1977         9.53     NaN
## YearF1978        11.58     NaN
## YearF1979        11.61     NaN
## YearF1980        12.24     NaN
## YearF1981        13.35     NaN
## YearF1982        11.87     NaN
## YearF1983        12.14     NaN
## YearF1984        12.45     NaN
## YearF1985        11.03     NaN
## YearF1986        11.23     NaN
## YearF1987        10.88     NaN
## YearF1988        11.42     NaN
## YearF1989        12.21     NaN
## YearF1990        10.85     NaN
## YearF1991        10.92     NaN
## YearF1992        10.67     NaN
## YearF1993        11.15     NaN
## YearF1994        10.88     NaN
## YearF1995        10.50     NaN
## YearF1996        10.38     NaN
## YearF1997        10.63     NaN
## YearF1998        10.63     NaN
## YearF1999        10.60     NaN
## YearF2000        10.32     NaN
## YearF2001        10.08     NaN
## YearF2002         9.61     NaN
## YearF2003         9.65     NaN
## YearF2004         9.71     NaN
## YearF2005         9.46     NaN
## YearF2006         9.57     NaN
## YearF2007         9.93     NaN
## YearF2008        10.16     NaN
## YearF2009         9.47     NaN
## YearF2010         9.92     NaN
## YearF2011        10.17     NaN
## YearF2012         9.84     NaN
## YearF2013         9.72     NaN
## YearF2014         9.80     NaN
## YearF2015         9.81     NaN
## YearF2016         9.93     NaN
## YearF2017         9.59     NaN
## YearF2018         9.16     NaN
## YearF2019         8.59     NaN
## YearF2020         8.47     NaN
## YearF2021         8.32     NaN
## GearCatBT         0.25     NaN
## GearCatGOV_CL     1.54     NaN
## GearCatGOV_GG     0.97     NaN
## GearCatOTT       -3.26     NaN
## GearCatPEL       -2.05     NaN
## GearCatTV         0.89     NaN
## 
## Smooth terms:
##            Std. Dev.
## sds(Month)      0.43
## 
## Dispersion parameter: 665.08
## Tweedie p: 1.61
## Matern range: 1.35
## Spatial SD: 0.37
## Spatiotemporal SD: 0.26
## ML criterion at convergence: 248325.051
## 
## See ?tidy.sdmTMB to extract these values as a data frame.
#g <- visreg::visreg(fit, xvar = "GearCat", gg = TRUE)
#plot(g)
p1=visreg::visreg(m1, xvar = "GearCat", gg = TRUE)
p2=visreg::visreg(m1, xvar = "YearF",gg = TRUE)
p3=visreg::visreg(m1, xvar = "Month",gg = TRUE)
p1

p2

p3

Predict onto the grid

grid=grid %>% filter(!Year %in% c(1967,1968))
grid$Year=as.integer(grid$Year)
grid$YearF=as.factor(grid$Year)
dat %>% group_by(GearCat) %>% summarise(n=n())
## # A tibble: 7 × 2
##   GearCat     n
##   <fct>   <int>
## 1 BAK_GG    520
## 2 BT      33421
## 3 GOV_CL  37009
## 4 GOV_GG   6353
## 5 OTT        87
## 6 PEL        88
## 7 TV      10896
grid$GearCat = 'GOV_CL'
# calculate the mode of month
# getmode=function(x){uniq=unique(x)
# uniq[which.max(tabulate(match(x,uniq)))]}
# getmode(dat$Month) Month set to 9
grid$Month = 9

p <- predict(m1, newdata = grid)

Spatial random fields

plot_map(p, "omega_s") +
  scale_fill_gradient2() +
  ggtitle("Spatial random effects only")
## Warning: Removed 4817 rows containing missing values (geom_raster).

Spatio-temporal random fields

plot_map(filter(p,Year%%5==1), "epsilon_st") +
  scale_fill_gradient2() +
  facet_wrap(~Year,ncol = 2) +
  ggtitle("Spatiotemporal random effects only")
## Warning: Removed 1169 rows containing missing values (geom_raster).

plot_map(filter(p,Year%%5==1), "exp(est)/1000") +
  scale_fill_viridis_c(
    trans = "pseudo_log",
    # trim extreme high values to make spatial variation more visible
    na.value = "yellow", limits = c(0, quantile(exp(p$est)/1000, 0.995)),'cod adult biomass \n (kg km^-2)'
  ) +
  ggtitle("Prediction (fixed effects + all random effects)",
    subtitle = paste("maximum estimated biomass density (whole time series) =", round(max(exp(p$est))/1000))
  )+facet_wrap(~Year,ncol=2)
## Warning: Removed 1169 rows containing missing values (geom_raster).